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SUMMARY 


An on-line technique is presented for the identification of rotor blade 
modal damping and frequency from rotorcraft random response test data. The 
identification technique is based upon a recursive maximum likelihood (RML) 
algorithm, which is demonstrated to have excellent convergence characteristics 
in the presence of random measurement noise and random excitation. The RML 
technique requires virtually no user interaction, provides accurate confidence 
bands on the parameter estimates, and can be used for continuous monitoring of 
modal damping and frequency during wind tunnel or flight testing. 

Results are presented from simulation random response data which quantify 
the identified parameter convergence behavior for various levels of random excita- 
tion. The data length required for acceptable parameter accuracy is shown to 
depend upon the amplitude of random response and the modal damping level. Random 
response amplitudes of 1.25 degrees to .05 degrees are investigated. The RML 
technique was applied to hingeless rotor test data from the NASA Langley Research 
Center Helicopter Hover Facility. The inplane lag regressing mode damping was 
correctly identified at different rotor speeds. The identification from the test 
data agreed with the simulation results and with other available estimates of 
frequency and damping. 


INTRODUCTION 

The determination of rotor blade modal damping for lightly damped rotors 
during wind tunnel or flight testing is an important safety requirement. Analytic 
methods to predict stability are often employed before testing, however the 
accuracy of analytic methods are not always reliable, particularly for new rotor 
designs or configuration modifications. Estimation of rotor blade mode damping 
from test data is thus required during wind tunnel testing where the possibility 
of low damping levels may occur. Current methods for estimation of modal damping 
from test data are based upon use of the Fast Fourier Transform (FFT) . This 
technique is often referred to as the moving block FFT method. This method has 
been applied to rotor data obtained either from an initial condition decay 
response obtained after removal of a forced excitation, or from random response 
data (ref. 1, ref. 2 and ref. 3). The random response data is first transformed 
to an initial condition decay response by calculating an ensemble average of 
segments of the random response time history. Application of the moving block 
FFT approach is then performed. This technique has been referred to as the moving- 
block/ randomdec method (ref. 1). 

During wind tunnel testing, it would be desirable if damping and frequency 
could be computed without the use of a forced input excitation and be computed 
in real time along with confidence bands. Although the FFT technique has 
received widespread acceptance by the rotorcraft community, practical usage of 
this technique requires considerable user interaction and is not a real time method 
Estimation of damping and frequency requires several minutes and continued user 



interaction. In addition, wind tunnel applications are usually done with a forced 
input excitation to obtain the decay response time history. 

The 'RMl, technique presented in this report requires virtually no user interac- 
tion, and provides on-line continuous estimates of damping and frequency from 
random response data. In addition, accurate confidence bands (+ la bands) are 
also recursively computed which indicate the required data length for accurate 
parameter estimation. The technique can be used on very small amplitude random 
response data which often is the case in wind tunnel tests. During wind tunnel 
testing it is common to have data continuously generated while steady trim condi- 
tions are achieved, and this data is suitable for parameter identification. 

Thus, the identification procedure need not interfere with other wind tunnel test 
objectives. Although the RML technique presented is developed for any number of modes, 
including closely spaced or overlapping modes, the simulation and test results 
presented in this report focus on the single mode case. 

Helicopters can encounter many different types of rotor dynamic instability 
during wind tunnel testing. Instability can arise from pure mechanical coupling 
as in ground resonance or from aerodynamic coupling in forward flight. Rotor 
blades are designed such that the rotor blade mode frequencies lie between 
multiples of the rotor normal operating RPM. This fact places the modes of 
Interest between a known band of frequencies and allows isolation of these modes 
with the use of a bandpass filter. (This is in contrast to fixed wing flutter 
where the modal frequencies are not as well defined.) Prior knowledge of the 
approximate frequency band allows for pretest selection of the bandpass 
filter passband, which in turn permits automated use of the identification 
algorithm. The identification method identifies both damping and frequency of 
the modes within the preselected passband. 

Recursive maximum likelihood methods for the identification of parameters 
from input/output data have been extensively researched (ref. 4 and ref. 5). 

However, the application to real-time identification of rotor blade damping has 
not been previously done. These real-time methods are based upon an input /output 
discrete transfer relationship or autoregressive moving average (ARMA) representa- 
tion. This representation leads to a simple, accurate, and numerically efficient 
identification algorithm. Other methods for real-time identification of para- 
meters such as least squares or instrumental variable produce either biased 
estimates in the presence of measurement noise or require extensive user pre- 
tuning for proper convergence. The recursive maximum likelihood procedure selected 
for the research reported herein has excellent convergence properties and requires 
no user pretuning. This technique is referred to as the RML2 method in reference 
4 and reference 5. 

Flutter boundary prediction has been performed in reference 6 for a fixed 
wing aircraft using nonrecursive maximum likelihood identification in conjuction 
with Jury's determinate method to predict the stability boundary. The present 
research differs from the techniques used in reference 6 in that recursive on- 
line maximum likelihood procedures are used. In addition, the objective in this 
research is to provide continuous on-line estimation of modal damping as wind 
tunnel testing is performed. 

The real time Identification procedure is described in this report. An 
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investigation is conducted using simulation data to determine the convergence 
characteristics of the algorithm and identification accuracy with random input 
excitation and measurement noise. The effect of bandpass filter bandwidth on 
identification accuracy is quantified. A lightly damped blade mode with 2% 
critical damping and moderately damped mode of 10% critical damping is used in 
the investigation. The algorithm is applied to hingeless rotor random response 
test data obtained from the NASA Langley Research Center Hover Test Facility. 
Test data is analyzed using the algorithm and damping and frequency estimates 
are compared with other available predictions. 
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SYMBOLS 


coefficients of the characteristic equation for a dynamic system of 
order n, j = 1, n, also defined in equation 2.9 
stability matrix appearing in linear state equation 2.5 

polynomial in the shift operator z”^, A(z) = 1 + aj^z”^ + a 2 z”^ + + ^n^'^ 

control matrix appearing in linear state equation 2.5 

coefficients associated with the control input in the ARMA representation 

polynomial in the shift operator z"^, B(z) = biz“^ + b 2 Z~^ H + b^z“^ 

cyclic control input to rotor, deg. 

filter bandwidth used in the digitial filter, Hz 

coefficients associated with the innovations in the ARMA representation 


polynomial on the shift operator z"^, C (z) = 1 + c-|Z ^ + C 2 Z ^ + + Cj^z ^ 


covariance of the parameter estimate 
covariance of the parameter estimate corrected for ndnwhite innovations 
natural logrithm basis, e = 2.71828 

generalized force scaled by generalized inertia for the ith mode 

control matrix appearing in the discrete state equation 2.7 

measurement matrix appearing in the linear state equation 2.7 

weights used in the digital filter 

discrete data point index, also mode number index 

complex variable j = also data sample shift index 

Kalman f liter gain 

gain matrix used in the recursive maximum likelihood algorithm 
negative log likelihood function 
likelihood function 

number of data points on either side of current point fitered by the 
digital filter 

vector of mode shapes at a radial station 
number of data points 


number of coefficients in ARMA model for a^^, bj_, c^; i = l,n 
parameter covariance matrix used in the recursive maximum likelihood 
algorithm 

initial parameter covariance 
nondimens ional rotor frequency 

generalized rotor coordinate associated with ith mode 
root mean square 
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Laplace transform variable s =» a 4- jo) 
time variable, sec, also discrete time index 
data sample time interval, sec 
control input 

discrete white gaussian random measurement noise 
loss function minimized in the RML algorithm 
discrete white gaussian random process noise 
zero-mean continuous white gaussian process noise 
state vector 

vector of sensitivity equations used in the RML algorithm 
measurement variable (or vector) 
digital filtered data 

z-transform shift variable, also modal response variable 

real and Imaginary parts of the complex variable Z,Z=Zp^+jZj 

innovations 

process noise gain matrix 
innovation covariance 
convergence parameter 
innovations 

discrete state transition matrix 
parameter vector 
initial parameter vector 

standard deviation, also used as real part of eigenvalue 

standard deviation of process noise and measurement noise respectively 

standard deviation of damping parameter estimate 

frequency, rad/sec. 

undamped natural frequency, rad/sec 

damped frequency associated with inplane mode, rad/sec 

cutoff and termination frequency used in the digital filter, rad/sec 

rotor speed, revolutions per min 

modal damping, % critical 

denotes transpose 

denotes time derivative 

denotes estimated quantity 

partial derivative with respect to 0 
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THEORETICAL BACKGROUND 


This section presents the problem description and overall approach used in 
the identification of blade modal damping and frequency from test data* The 
real-time maximum likelihood identification algorithm is then presented followed 
by the bandpass filter used to separate the desired modes from the measured test 
data. 


Problem Description and Overall Approach 


The real time identification of helicopter blade damping and frequency is 
based upon the blade modal response to random excitation. This response is 
assumed to be modeled by a set of coupled second order equations with a single 
output measurement occurring at discrete time increments. The rotor blade response 
contains frequencies occurring at multiples of the rotor frequency (Ip, 2p, etc.) 
and the various blade modes. Bandpass filtering is employed to separate the main 
modes of interest from the rotor Ip, 2p, etc. frequencies and other modes. The 
resulting response data can then be considered to be a low order model of N 
coupled second order equations. Rotor blade response in this passband will 
typically consist of from one to three blade modes. 


The general response of the ith blade mode of the rotor in the rotating 
system can be represented as 






% 


(JL) 


q. = F. 


( 2 . 1 ) 


where q^ is the generalized coordinate and is the generalized force (scaled 
by the generalized inertia) of the ith mode. The vertical deflection, Z(t) at a 
radial station is the sum of the modes 


// modes 

z(t) = E qi 

i=i 


( 2 . 2 ) 


where is a vector of mode shapes at station r. The variable Z(t) can also be 
used to represent the blade bending moment at station r as measured by strain 
gauage instrumentation. Defining qj^ = ^i^i substituting into Eq. (2.1) and 

Eq. (2.2) eliminates the mode shape vector from the vertical deflection Z. 
Equations (1.1) and (2.2) thus becomes 

ii + ^ = h (2.3) 

and 

// modes _ 

z(t) = E q^ ( 2 . 4 ) 

i=l 

where is scaled by the corresponding mode shape 

A canonical representation for the deterministic single input-single output 
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system described by Eq. (2.3) and Eq. (2.4) (ref. 7) can be described in state 
form as 


X = Ax + B'u 


(2.5) 


where 



( 2 . 6 ) 



B' = 


^2 


K-1 


b' , 
n-2 


and X is a n x 1 state vector of normalized coordinates and velocities, u is 
the single input and y is the single output. 

In the presence of random process and measurement noise the continuous system 
Eqs. (2.5) and (2.6) can be represented by the discrete stochastic system of linear 
difference equations 


x(i+l) = (j)x(l) + Gu(i) + rw(i) 
y(i) = Hx(i) + v(i) 


(2.7) 


where i is the index for the discrete sample time. The discrete time process 
noise w(i) and measurement noise v(i) are assumed to be zero mean white sequences 
with covariance matrices Q and R respectively. 

By use of the Kalman Filtering theorem (ref. 8 and ref. 9) the stochastic 
model of Eq. (2.7) can be represented by an equivalent form as 
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x(i+l) = (p x(i) + Gu(l) + Kv(i) 
y(i) = H x(i) + v(i) 


( 2 . 8 ) 


where x(i) denotes the conditional mean of x(i) given measurements y(i-l), y(i-2), 

« . • , and {v(i) } is a sequence of independent equally distributed random variables 
with zero mean and covariance The representation of Eq. (2.8) is an innova- 

tion representation (ref. 8 and ref. 9) of the process and the sequence {v(t)} 
is called the Innovations. In steady state the filter gain matrix K is constant 
and combined with the innovations covariance X^, replaces the noise statistics Q 
and R of Eq. (2.7). The innovations sequence v(i) replaces both the process noise 
and measurement noise. 

The innovation model of Eq. (2.8) takes on a particularly simple form for a 
single input-single output (SISO) model. The discrete SISO form or autoregressive 
moving average (ARMA) representation is obtained using the shift operator q“"^, where 
q*“lx(i+l) = x(i). Application of the shift operator to Eq. (2.8) yields the SISO 
ARMA representation 

n-1 n 

2 a.y(i-j) = V {b.u(j-l) + c.V(i-j)} (2.9) 

j=0 ^ ^ ^ 


where a^^l and Cq = 1, n is the order of the autoregressive part, and also the 
order of the moving average part. In this description there is no concept of 
state and a n state vector produces n lags y(i-j), j=0, 1, 2, ...n-1. 


The ARMA representation of Eq. (2.9) provides a complete description of the 
SISO stochastic process of Eq. (2.7). The unknown parameters of the model in 
Eq. (2.9) are the parameter a j , b j , Cj and the scalar innovation covariance. The 
unknown parameter vector 0 is 



( 2 . 10 ) 


along with the scalar innovation covariance denoted by X^. The left-hand side of 
Eq. (2.9) contains the coefficients of the characteristic equation. 


The damping and frequency of the n-modes are obtained from the roots of the 
characteristic equation. 

Z^~^ + = 0 ( 2 . 11 ) 

where Z is the discrete Z transform variable. 

The general ARMA model given by Eq. (2.9) can thus be used to identify the 
coefficients a^ , j=l, n-1 of a general n^^ order system from the single measure- 
ment y. The stability of the system is determined from the roots of the discrete 
characteristic equation. To obtain the continuous damping and frequency the 
mapping from the discrete z-transform domain to the continuous s-transform Laplace 
transform domain is then required. 
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The ARMA representation of Eq. (2,9) is more suitable for on-line identifica- 
tion than the complete state space model of Eq. (2.7) or the state innovations 
model of Eq. (2.8). The ARMA representation has fewer parameters and the noise 
statistics are completely characterized in the moving average part of Eq. (2.9) 

(i.e. the c coefficients) and innovation covariance Thus in order to identify 

the coefficients of the characteristic equation (i.e. the a coefficients) , it is 
also required to identify the c coefficients and innovation covariance This 

then provides a complete description of the stochastic process, from which the 
modal damping and frequency can be determined. Also, if a known control input is 
applied the b coefficients of Eq. (2.9) must also be identified. 

The overall identification procedure is shown in block diagram form in figure 
1. The data y (and control input u(i) if available) is discretized and bandpass 
filtered to obtain the desired rotor blade modes. The filtered data is passed 
through the recursive maximum likelihood identification algorithm (discussed in 
the next section) and the identified parameters represent the coefficients of 
the discrete characteristic equation. Discrete roots (eigenvalues) are computed 
from the characteristic equation. The modal damping and frequency are then obtained 
by transforming the discrete roots to the continuous Laplace transform domain. The 
RML algorithm computes both the parameters and standard deviations of the parameter 
estimates as noted in figure 1. 

The discrete roots are obtained from the characteristic equation by use of any 
standard root finding algorithm of a polynomial equation. Transformation to the 
continuous domain is accomplished by defining the discrete complex root Z and 
continuous root s 


z = Zr + j 


( 2 . 12 ) 


s = a + j 0 ) 


and using the discrete to continuous transformation 


sT 

e 


Z 


(2.13) 


(2.14) 


where T represents the discrete sample time interval. 


Substituting Eq. (2.12) and (2.13) into Eq. (2.14) yields 


^0T + jwT 


\ + j 


(2.15) 


Expressing Eq. (2.15) 
hand sides yields 


in terms of amplitude and phase and equating right and left 


-4 + 4 

-1 ^1 

ojT. = tan 


(2.16) 

(2.17) 


Taking the natural logrithm of Eq. (2.16) and solving for 0 and OJ in Eq. (2.16) and 
(2.17) yields the continuous damping and frequency 
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(2.18) 


0 - -ly In (z2 + 

2 

u=itan~^ (2.19) 

The undamped natural frequence and damping ratio ^ are related by 

= a (2.20) 

( 2 . 21 ) 

can be computed with the use of 

( 2 . 22 ) 


w = ojn (1-C^)^^^ 


The damping ratio and undamped natural frequency 
Eq. (2.18) through Eq. (2.21) 

1/2 


/ 2 2. 
0 ) = (o) - a ) 
n 



(2.23) 


If there is only one mode, the root finding algorithm is not required and the 
damping and frequency can be obtained directly from the a coefficients a^ and a2 
using 

a = In (1/32) (2.24) 


0 ) = -| cos ^(-ai/2/a2) (2.25) 

The next section describes the details of the recursive maximum likelihood 
identification algorithm for computing the coefficients of the ARMA represention 
of Eq. (2. 9) . 

Recursive Maximum Likelihood Identification Algorithm 

The origin of this method is based upon the off-line maximum likelihood method 
of Astrom-Bohlin (ref. 10). The recursive version is described in Soderstrom (ref. 
11). This method is selected for the rotorcraft blade mode identification due to 
its excellent convergence properties. Reference 4 and Reference 5 describes the 
recursive version in comparison to other methods. 


The off-line maximum likelihood version attempts to minimize the negative log 
likelihood function 


L’ = -log L(0,X) 


12 

— y y) e + N/2 log A + N/2 log 2 tt 
2A ^1 ^ 


where 


0 — ^2* * 


.. a^, b^, 


n 


^1’ *^2’ 


(2.26) 


(2.27) 
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and represents the innovation covariance. The ARMA model of Eq. (2.9) is 
rewritten using the operator C(z) , A(z) and B(z) 


C(z) = A(z) - B(z) , t = 1,2,. ..N (2.28) 

where the innovations is denoted by and t denotes the discrete sample time. 

Since the prediction error Eq. (2.28) is nonlinear in the parameters c, a 
iterative solution is required to minimize the likelihood function L*. 

2 

Optimization of L’ with respect to A yields 

9^9 

X = 1/N X) (2.29) 

t=i 

Therefore maximum likelihood identification of 0 is equivalent to minimization 
of 

N 

V(9) = E (2.30) 

t=l 

subject to the prediction error model 

= [C(z)]“^ [A(z) - B(z) u^] (2.31) 

An approximate solution is obtained by linearization of 


St = e(0) + || (0) (0 - 0) 


(2.32) 


where the linearization is performed about the estimate 0. 

Eq. (2.32) is a linear prediction error equation and therefore the minimization 
of Eq. (2.30) subject to Eq. (2.32) can be obtained using the recursive least square 
solution (ref. 12) 

Therefore the recursive solution is 


+ "t+1 Vl 

(2.33) 

*'t+l ■ ''t + ?+l ''t 

(2.34) 

Vi - -\+i ?+i 

(2.35) 


where is computed from 


10 
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(2.36) 


where (t) = e(t). 

The sensitivity =• ^ is computed according to 



(2.37) 


The factor in Eq. (2.34) and (2.35) is used to improve convergence rate 
(ref. 4) and is obtained from 

= C X*^ + 1 - C (2.38) 

c o c o 

where X° = .95 and C = . 99 . The values of X° and C are the same as those used in 
reference 4 which were found to yield good convergence. 
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To track slowly time varying parameters the convergence factor is set to 
a constant value less than unity for all time. This is equivalent to exponential 
weighting of the least square criterion Eq. (2.30) and is fully discussed in 
reference 4. 

Eqs. (2.33) through Eq. (2.38) constitutes the recursive maximum likelihood 
algorithm. 

Bandpass Filter 

The selection of the desired modes for identification is accomplished with a 
zero phase shift digital filter, due to Graham (ref. 13). This filter is not 
recursive but uses numerical convolution of the raw data to obtain the filtered 
data. This filter can be used in near-real time by shifting the filter weight 
function hj one data point ahead as each new raw data point is received. The 
actual filtered output requires a finite segment of data ahead in time of the 
current data point being filtered. This results in a slight time delay 
before the estimated coefficients are available. Using a sample rate of .01 
seconds this represents a delay of between 1 and 2 seconds before the filtered 
output is obtained. A recursive bandpass filter such as a n^^ order Butterworth 
filter could be used to eliminate this delay, however the Graham filter was 
selected for the results presented here due to its excellent filter properties 
and availability. The 1 to 2 second delay was not considered significant for 
the real time monitoring of blade modes. 

The filtered output y^ at data point i obtained using the numerical 
convolution is 



h. y... 

J J+i 


(2.39) 


where yj+i is the raw data which requires M data points ahead and M data points 
behind point i. The weight function hj is computed from 


h. 

3 


IT 

2jT 


(sin ca^jT + sin o) jT) 

7T^ - 


3 


-M, . . . ,+M 


(2.40) 


j ^ 0 


and 


h. 

3 


0)^ + to 
t c 

27T 


3 = 0 


The cutoff frequency and termination frequency 0)^. are selected to achieve the 
desired passband about a center frequency of zero (low pass filter) . The weights 
h^ are then normalized by the sum of all the weights. To achieve bandpass filtering 
about a center frequency U)q, the weights are shifted by multiplication of h^ by 
cos ^9* (2.39) and (2.40) describe the low pass filter (ref. 13) and^ the bandpass 

filter is obtained by multiplication of the weights hj by cos cOqJT and is described 
in reference 14. 
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The use of the bandpass filter on both simulation data and test data is 
presented in the next section along with the recursive maximum likelihood iden- 
tification results. 
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RESULTS FROM SIMULATION ANALYSIS 

Before analysis of test data is performed it is important to assess the 
characteristics and limitations of the identification method by simulation. Of 
particular importance is the convergence characteristics. (i.e. the quantifica- 
tion of parameter accuracy as the data is processed) The amplitude of the 
random response, the initial parameter start-up, and the digital filter are 
investigated as to their effect on identified parameter accuracy. In addition, 
the standard deviation of the parameter estimate as computed in the algorithm 
must be evaluated as to its validity as a measure of identified parameter 
accuracy. 


Simulation Model 

The model used in the simulation analysis is intended to represent two situa- 
tions. First, a moderately damped mode (C = 10% critical damping) representative 
of an articulated rotor or highly damped hingless rotor is investigated. Second, 
a lightly damped mode (? = 2% critical damping) representation of a hingless rotor 
is investigated. In both cases, the simulation model represents a single uncoupled 
mode. Such a mode can occur as an isolated blade mode or a lightly damped "system” 
mode such as the in-plane regressing mode typical of ground resonance. Although 
the identification procedure is applicable to any number of coupled modes, the 
results presented are for the single mode case which is perhaps the most important 
for rotor stability assessment. In addition, the modal identification procedure 
can be repeated any number of times on the same data, in each case a different 
mode isolated with the bandpass filter. 

The model used in the simulation is a state space second order system and 
is represented by 


Xi = X2 


= -0) X- - 2^0) + w(t) 

z n 1 n z 


(3.1) 

(3.2) 


where x^ is the blade modal response and w(t) is zero-mean white gaussian process 
noise. 

The modal measurement is obtained from 

y = x^ + r (3.3) 


where, r is zero-mean white gaussian measurement noise. 


The 

discrete 


stochastic model given by Eq. (3.1) through Eq. (3.3) is simulated in 
form by 


x^(i+l) = Xj^(i) + X 2 (i) At (3.4) 

X2(i+1) = At Xj^(i) + (l-2?w^ At) X2(i) + w^(i) 


14 



and 


y(i) = x^(i) + r(i) 


(3.5) 


where is zero-mean white discrete process noise’ and r is zero-mean white discrete 
measurement noise. 

In order to obtain acceptable integration accuracy using the discrete model 
of Eq. (3-4) and Eq. (3.5), the integration step size At should be at least 10 
times smaller than the data sample time between the i and i+1 sample. For all 
simulation runs At was selected as At = .002 seconds and the measurement y(i) was 
made every At = *02 seconds. 

The maximum likelihood algorithm identifies the coefficients of the equivalent 
ARMA representation as discussed in the previous section. The ARMA representation 
used is 




(3.6) 


As discussed previously, Eq. (3.6) is a single input-single output discrete 
model and is equivalent to an innovations or Kalman filter representation of the 
system described by Eq. (3.4) through Eq. (3.5). The filter gains are embedded 
in the c coefficients of Eq. (3.6). The a coefficients are the coefficients of 
the characteristic equation from which the damping, and frequency, can be 

computed as shown in Eq. (2.20) through Eq. (2.25) of the previous section. 

Moderate Damping Case, ^ = 10%, = 23.15 rad/sec 

The simulation model shown by Eq. (3.4) and Eq. (3.5) was used to generate 
5000 data samples. The measurement noise standard deviation used is = .1 and 
various levels of process noise were investigated. The data used in the identifi- 
cation algorithm is unfiltered. The initial parameter start-up is 

0(0) = {-1.5 .80 0. 0.} (3.7) 


The parameter a2 = *80 represents a critical damping of greater than 20%. 

The initial parameter covariance used is 

P(0) = {.1 .1 1.0 1.0} (3.8) 

Three different levels of process noise were simulated (a^ = 1.0, .2, and .04). 
This corresponds to a state response with RMS levels of 1-25 deg, .25 deg, and .05 
deg, respectively. Identified parameter convergence is shown in figure 2 through 
figure 4 . f or the three levels of random excitation. As shown in the figures, con- 
vergence to the true parameter values occurrs in each case, with slower convergence 
achieved for the smaller levels of process noise excitation. The ARMA coefficients 
are shown on the left hand scale and the damping ratio is shown on the right hand 
scale of the plot for coefficient a 2 * The damping ratio is computed according to 

n 


The identified c coefficients account for 
If the ratio of process to measurement noise is 


the presence of measurement 
large as in figure 2 (cr^/o^ 


noise. 

= 100 ) , 
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2 2 

then the c coefficients are small- As the ratio a /a decreases, the c coefficients 

^ ^ 2 2 

increase in magnitude as shown in figure 4. In the limit as approaches then 

c approaches zero, which reduces to a least square estimator for the a coefficients. 
Thus, the c coefficients in the maximum likelihood method accounts for the presence 
of both measurement and process noise by automatically ^tuning" the Kalman filter in 
proportion to the ratio of process to measurement noise . 


Effect of Initial Parameter Estimate . The parameter identification results 
just presented for the 10% critical damping case is repeated for different initial 
parameter estimates. The second set for the initial parameter estimate and covariance 
is 


0 ( 0 ) = {0 0 0 0 } ( 3 . 10 ) 

P(0) = {1.0 1.0 1.0 1.0} (3.11) 


The identified parameter convergence is shown in figure 5 and figure 6 for 
process noise levels of = .2 and = .04, respectively. Again, identified 
parameter convergence is excellent. In all cases studied the final convergence 
was found to be independent of the initial parameter start-up. 

Lightly Damped Case, ^ = 2%, co =23.15 rad/sec 

n 

The simulation was repeated using a critical damping of ? = 2%. This case is 
representative of a. lightly damped hingeless rotor. Two levels of process noise 
were used (a^ = .2 and = .04). This random excitation produces a RMS random 
response of .75 deg and .15 deg, respectively. The identified parameter convergence 
for these two cases is shown in figure 7 and figure 8. Parameter convergence is 
excellent and occurs more rapidly than for the 10% damping case (compare fig. 7 
with fig. 3 and fig. 8 with fig. 4). The RMS of the random response is larger for 
the low damping case, as is expected, and the damping is more rapidly identified* 

Data Length Required For Accurate Identification 

Based upon the limited number of simulation runs performed an estimate of the 
number of data samples required for accurate identification can be approximated. As 
was shown in the simulation study, the required data length for convergence depends 
on the process noise level and damping value. Since the actual process noise level 
is not known in a wind tunnel test program, the RMS amplitude of the random response 
is presented rather than process noise level (o^) , Figure 9 presents a graphical 
summary of the simulations performed and shows the number of data samples required 
to identify damping ratio (% critical) to within 10% error. This is plotted vs the 
RMS amplitude of random response in degrees. Two different damping ratios (C = 10% 
and 5 = 2%) are presented. The results shown in figure 9 are obtained for a measure- 
ment noise level of Oj- = .1 degrees, O)^ = 23.15 rad/sec, and At = .02 seconds sample 
time interval. 

The number of samples required is very dependent upon the RMS amplitude and 
the damping ratio. Typically, the number of samples required ranges from 1000 
samples for low damping (2%) and large RMS amplitude response (1.25 deg) to greater 
than 5000 samples for high damping (10%) and low RMS response (.05 degrees). The 
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results in figure 9 should be used only as a guide since the results are presented 
for one value of measurement noise (a^ = .1), using unfiltered data, and a sampling 
rate of At = .02 seconds. 

The results of figure 9 show that identified parameter convergence requires 
between 20 seconds to 100 seconds of data (At = .02) depending upon the amplitude 
of random response and system damping level. A significant finding is that damping 
can be accurately identified from random response levels as low as .05 degrees in 
the presence of .1 degrees measurement noise. The smaller the RMS response, the 
more data samples required for parameter convergence. 

Digital Filter Selection 

A bandpass filter is required to remove unwanted modes and 1 per rev contamina- 
tions from the data before identification is performed. This section demonstrates 
the effect of the bandpass filter on identification accuracy for various filter 
designs. Figure 10 shows the amplitude vs. frequency characteristics of four 
digital filters used to separate the desired mode. The mode of interest, o)^, occurs 
at 3.68 Hz (23.15 rad/ sec) and the 1 per rev (Ip) signal is shown at 11.17 Hz. 

Filter number 1, 2 and 4 of figure 10 are bandpass filters, whereas filter number 3 
is a low pass filter. 

The simulation data was generated with a process noise level of - .2 for 
the 2% critical damping case. Four separate runs were made and identified para- 
meter convergence was assessed. Figure 11 shows identified parameter convergence 
for filter number 4. The identified damping value is somewhat biased from the 
true value of ? = 2%. Filter number 1, 2 and 3 showed a larger bias, indicating 
that some of the signal (i.e. the random response) was removed by the filter. The 
smaller the filter bandwidth, the more signal removed and thus, the larger the 
error. Figure 12 shows the effect of filter bandwidth on identified damping ratio 
for the four filters used. 

It is shown that a filter bandwidth between 8 and 10 Hz is required to achieve 
acceptable accuracy. Normalizing the frequency by the natural frequency, the following 
approximate frequency range should be selected for the filter assuming a white noise 
random excitation 

.2703^ 5 0) £ 2.4co^ (3.12) 

The relationship in Eq. (3.12) is approximate and in practice the filter band- 
width should be chosen as large as possible within the 0 to Ip frequency band for 
modes below 1 p. For modes above Ip but below 2p the filter should be selected 
between Ip and 2p frequencies. 

For a mode which lies very close to the Ip frequency, it may be required to 
remove the Ip frequency with a separate estimator before using the bandpass 
filter. 

Figure 11 shows that the identified a coefficients are somewhat biased by 
the filter. The c coefficients are significantly biased. This is because the 
c coefficients are influenced by the ratio of process noise to measurement noise 
which assumes the random noise is white noise. Since filtering produces colored 
innovations, this is reflected in the identified c coefficients. 
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Covariance of the Parameter Estimates 


Real time identification of modal damping and frequency is not complete with- 
out knowledge of the uncertainty in the parameter estimate. The recursive maximum 
likelihood algorithm provdies a measure of the parameter covariance by use of the 
Cramer-Rao lower bound. A discussion of theorems relating the accuracy of maximum 
likelihood estimate to the Cramer-Rao lower bound is given in reference 12. The 
Cramer-Rao bound is a lower bound on estimation accuracy which is achieved by the 
maximiim likelihood estimator asymptotically. An estimate of the parameter covariance 
is obtained by 

cov(0jj) = Pjj xj (3.13) 

where P^ is obtained from Eq. (2.35) and the innovations covariance is computed 
from 



Since, in practice, the noise characteristics of the innovations is often 
found to be "colored" rather than white, the parameter covariance as given by Eq. 
(3.13) must be corrected for non-white innovations (ref. 15). This can be 
accomplished using 

^(0) = cov(0) • 

where, cov(0) is the covariance corrected for non-white innovations and BW is the 
bandwidth of the colored noise which is set equal to the digital filter bandwidth 
of the bandpass filter used on the data. The correction term is approximate and 
is found to work well in practice as discussed in reference 15. 

Since cov (0) represents the covariance of the coefficients of the charac- 
teristic equation, a transformation is required to obtain the covariance of the 
damping. Since the damping ratio is related to the coefficient a 2 by 


^ -2?(%At (3.16) 

a^ = e *■* 

the covariance of $ can be obtained by linearization of Eq. (3.16) to first order 
as 


£„ = 1 - 2CW At 
z n 

The covariance of C can then be computed from Eq. (3.17) or 

cov(8) = cov(a„)/(2 At w )^ 
z n 

The standard deviation of ? is thus 


(3.17) 


(3.18) 


v'cov(S-) 

ap = — (3.19) 

^ 2 At w 

n 
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Eq. (3.15) is used to compute the standard deviation of the parameters in 
the ARMA model (e.g. the a- coefficients) and Eq. (3.19) is used to estimate the 
standard deviation of the estimated damping. Figure 13 shows the identified para- 
meter convergence and standard deviation of the estimate for and a 2 for a 
simulation damping ratio of 2% and .2. The estimated standard deviation pro- 

vides an excellent measure of uncertainties for the identified parameters. 

Figure 14 shows identified parameter convergence and + 2a bands for a simula- 
tion damping of 2% with process noise = .04. The la bands were slightly 
optimistic while the 2a bands always contains the true damping value. Conver- 
gence is slower for figure 14 than figure 13 due to the smaller RMS of the random 
response (a^ = .04 vs a^ = .2). In all simulation studies performed the + la 
band or + 2a band provided an accurate measure of the parameter uncertainty. 
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RESULTS FROM TEST DATA 


Test Data Description 

A limited amount of test data obtained from a 4-bladed hingeless rotor test 
in the NASA Langley Research Center (LRC) Helicopter Hover Facility in 1981 was 
furnished by Mr. William T. Yeager, Jr. and Mr. Wayne R. Mantay. This data is 
used to identify in-plane mode damping and frequency. The rotor was operated in 
hover on the Aeroelastic Rotor Experimental System (ARES) over a rotor speed range 
of 0 to 618 RPM. The Helicopter Hover Facility is used for hover testing on the 
ARES before testing in the LRC Transonic Dynamics Wind Tunnel. Two RPM condi- 
tions, 250 RPM and 618 RPM, were investigated. 

Of particular concern during rotor speed increase is the potential for ground 
resonance which occurs near the rotor in-plane regressing mode frequency (f2-oa^) . 
Since the in-plane lead-lag angle measurement occurs in the rotating system the 
bandpass filter was chosen to be centered around the rotating in-plane frequency 
ca^. The 250RPM condition was observed to have a large resonance close to the 
rotor RPM, excited purely .by random turbulence. The in-plane magnitude response 
was approximately + 1 degree. Because of this large resonance the damping was 
expected to be quite small and a forced control input excitation would be dangerous 
to apply. The recursive maximum likelihood approach to identify modal damping from 
random response is thus ideally suited for this RPM condition. 

The 618 RPM condition was considerably more damped and the modal Response 
due to inherent random excitation was less than .05 degrees. It is expected 
that a longer data record length would be required for the 618 RPM condition 
for accurate identification because of the small amplitude response. 

The test data used for the 250 RPM and 618 RPM conditions are shown in figure 
15. Figure 15 shows only a portion of the 50 seconds of test data record length 
used in the identification- The lead-lag measurement for the 250 RPM case of 
figure 15 shows the + 1 degree resonance response occurring at in-plane frequency 
0)^. The sinusoidal response shown in figure 15 (bottom figure) for the measure- 
ment of lead-lag at the 618 RPM condition is due to the Ip frequency of the rotor 
and is not the in-plane mode. The large Ip response occurs at this condition due 
to a non-zero trim collective and cyclic control (0 q = 8°, - 4°). The in-plane 

mode amplitude is less than .05 degrees which will be discussed in more detail later. 

Identification Results at 250 RPM 

The data at the 250 RPM condition shown in figure 15 (top figure) was processed 
with bandpass filter number 4, the characteristics of which are shown in figure 10. 
During the process of transforming the analog data tape to digital format a portion 
of the data was distroyed. The region of bad data occurred between the 1200th and 
1600th samples. The complete data length processed consisted of 2500 samples at .02 
seconds per sample rate. 

The filtered data was used in the recursive maximum likelihood algorithm and 
the ARMA model coefficients identified. The results are shown in figure 16 for 
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the coefficients of the characteristic equation, aj^ and a.2 and the damping and 
frequency estimates are shown in figure 17. The region of bad data shown in 
these figures was also recursively processed by the algorithm. Identified 
parameter convergence is more accurate before the bad data region is encountered 
(i.e. up to 1200 samples). 

Convergence is shown to be excellent in figure 16 and 17 as indicated by the 
+ 2a accuracy bands shown in the figure before the region of bad data is encountered. 
The identified damping ratio at 1200 samples is 1.4% with 1 standard deviation equal 
+ .7%. The rapid convergence rate is a result of the large amplitude random re- 
sponse (+ 1 degree) as observed in the test data. The convergence rate shown in 
figure 16 and 17 is consistent with that expected from the simulation studies. 

For example, referring to figure 9, it is shown that perhaps 1500 data points are 
required to achieve the desired accuracy. Ideally, it would be desirable to pro- 
cess the full 2500 samples of data, however, because of the bad data region, the 
parameter estimate was shifted as shown in figure 16 and 17. Notice however, 
even with bad data, the algorithm is performing properly, in that the parameter 
estimates are attempting to recover as good data becomes available. 

The identified frequency is shown in figure 17 and the estimated value at 
1200 samples is 26.3 rad/sec. This value is nearly the same as determined by 
that obtained from shake test data and analytic predictions. 

The identified damping and frequency parameters shown in figure 17 demonstrates 
that the identification procedure is performing successfully, and is in excellent 
agreement to that obtained in the simulation studies. 

Identification Results at 618 RPM 

The hover test data obtained at the 618 RPM condition is shown again in 
figure 18 with an expanded scale. Figure 18 shows the unfiltered data (top 
figure) and the bandpass filtered data (lower figure) . The unfiltered data 
shows the dominant response is the Ip rotor frequency (+ .3 degrees). The 
filtered data removes the Ip frequency and retains the desired inplane mode 
which occurs at = 38.83 rad/sec. This mode is clearly shown in the lower 
figure of figure 18. The frequency of this mode appears to change frequency 
throughout the data length and is due to the random excitation driving the rotor 
system. The RMS amplitude of response is approximately .03 degrees as shown in 
the lower figure. Referring back to figure 10, the simulation study has shown 
that 5000 data points at a sample rate of .02 seconds (100 seconds of data) are 
required for acceptable accuracy for a RMS random response of .05 degrees. 

Since only 45 seconds of data are available at the 618 RPM condition and the RMS 
amplitude is -03 it is expected that convergence would require approximately twice 
as much data- This was in fact the case and the convergence plots are shown in 
figure 19. The identified damping ratio’ is shown in the lower figure and the 
+ 2a accuracy band shows a longer data record is required for convergence. The 
identified damping is estimated after 45 seconds of data (4500 samples at .01 
seconds sample rate) to be ^ = 9% with 1 standard deviation of + 3.5%. Indepen- 
dent estimates of the damping at this condition using the log decrement approach 
(visual inspection of an initial condition decay response was used) during the 
hover tests reveal a damping estimate of approximately 4%. The results of figure 
19 thus contain this value within the accuracy band (+ 2a) as shown in this figure. 
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The identified frequency convergence is also shown in figure 19 (top figure) 
and the estimated frequency is O) = 38.83 radians /sec. This value is close to 
that obtained from independent shake test data and analytic predictions. The 
frequency parameter is identified more rapidly than the damping and requires 
less data length. This was found to be the case in the simulation study per- 
formed prior to analyzing the test data. 

The identification results at the 618 RPM condition have correctly identified 
the frequency and approximately twice the data length available is required to 
identify the damping with acceptable accuracy. This conclusion was predicted by 
the simulation studies which show that 100 seconds of data is needed (twice the 
data length available) to obtain the required accuracy of the damping estimate. 

Table 1 provides a summary of the identified damping and frequency at both 
the 250 and 618 RPM conditions. Comparisons are made with other available pre- 
diction methods. Shake test data and analytic simulation is referred to in the 
table as Analytic Method. The Log Decrement method is based upon visual inspec- 
tion of the amplitude decay resulting from removal of a forced excitation. The 
damping values obtained from the Analytic method are estimates based upon the 
estimated rotor blade damping. The actual rotor blade damping value was not 
precisely known (1% critical damping was used in the simulation). The fixed 
system damping and frequency are for the in-plane regressing mode and the 
rotating system values are shown for the maximum likelihood method and log decre- 
ment method. The fixed system values are shown for the analytic method and 
rotating system values are obtained via transformation as noted in the table. 

From the data available at the 250 RPM and 618 RPM condition the recursive 
maximum likelihood approach performed as expected from the simulation studies. 

Identified damping convergence to acceptable accuracy depends upon the 
magnitude of the RMS random response and the level of modal damping. The larger 
the random response the faster convergence. Typically from 30 to 100 seconds of 
data are required for acceptable accuracy. For + 1 degree random response and 
low damping (1.4%) 30 to 50 seconds of data are required (250 RPM condition). 

For + .03 degree random response and damping approximately 4%, approximately 
100 seconds of data are required (618 RPM condition). 
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CONCLUSIONS 


A recursive maximum likelihood identification procedure has been presented 
for the estimation of helicopter rotor blade modal damping and frequency from 
random response data. The algorithm was shown to have excellent convergence 
properties, provide accurate estimates of parameter covariance, and require 
virtually no user interaction. As a result, the RML method can be used for 
continuous estimation of damping and frequency during wind tunnel (or flight) 
testing. 

Simulation studies have demonstrated that convergence to true parameter 
estimates occur even for very small amplitude RMS random response. RMS random 
response amplitude of .05 degrees to 1.25 degrees, which are typical of that 
expected in wind tunnel data, resulted in successful identification. Data record 
lengths required for accurate identification range from 20 seconds to 100 seconds 
of data depending upon RMS amplitude of random response and damping level. These 
data lengths are compatible with typical wind tunnel test data taken at steady 
rotor trim conditions. Since a forced input is not required, damping and frequency 
can be identified without interfering with other wind tunnel test objectives. The 
estimated parameter covariance was shown to provide accurate confidence bands on 
the damping, thus providing a indicator as to length of data required at a given 
test condition. Since identification parameter convergence was shown to be 
independent of initial parameter start-up values, the RML procedure can be used 
with virtually no user interaction. This proved to be the case on all data pro- 
cessed in this investigation. 

The hover facility test data identified parameter convergence results obtained 
were consistent with the simulation studies. At 250 RPM, identified in-plane mode 
parameters are t = 1.4% critical damping and = 26.3 rad/sec. The + la estimate 
on ^ is .7% using 24 seconds of data. The identified frequency is in agreement 
with independent shake test data results. Accurate estimates of damping from other 
independent tests were not available. At 618 RPM, the identified frequency was 
air = 38.83 rad/ sec which is in agreement with independent shake^test data. The 
identified damping is estimated using 45 seconds of data to be C “ 9% with 
+ la = 3.5%. The large standard deviation shows that a longer data record is 
required which is consistent with the simulation results which showed that twice 
the data length is needed to obtain an accurate estimate. For both the 250 RPM 
and 618 RPM the RML identified parameter convergence performed as predicted by 
the simulation studies. 
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RECOMMENDATIONS 


The emphasis of this investigation was on a single mode response separated 
from other modes with a bandpass filter. Although the techniques developed are 
applicable to any number of over-lapping modes, simulation studies on test data 
analysis have only been performed on a single mode case. Since rotorcraft rotor 
instability can result from coupled modes it is recommended that simulation 
studies analogous to that performed in this report be repeated for coupled mode 
response of two and possibly three modes. In addition, criteria should be 
investigated to automatically determine the required number of modes in the data. 

The RML technique contains provision to track time varying parameters. 

Simulations should be performed to quantify this capability. This would be 
useful for monitoring damping changes while going from one steady rotor trim 
condition to another trim condition. The RML technique contains provision for identi- 
fication of data generated by known forced control inputs and free response data. 

Both these situations have not been investigated via simulation. Simulation studies 
should be performed to quantify this capability. 

Finally, further test data applications should be performed. Both on-line and 
off-line applications would provide a data base for further assessment of the RML 
technique. 
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Table 1. - Test Data Identified Damping and Frequency and 
Comparisons with Analytic and Log-Decrement 
Predictions 


Estimation 

Rotating 

System 

Fixed 

System 

Method 


^R 

0) 

‘’F 


Rad/Sec 

% Critical 

Rad/Sec 

% Critical 

^ = 250 RPM 
Analytic Method 

26.2^^^ 

0.55%^^^ 

.8 

18% 

Log Decrement 

- Not Used^^^- 

— 

— 

Max. Likelihood 

26.3 

1.4%(la=.7%) 

, 

— 

0 , = 618 RPM 
Analytic Method 

37.63^^^ 

0.33%^'^^ 

! 

25.9 

.5% 

(2) 

Log Decrement 

— 

4% 

1 

1 

— 

Max. Likelihood 

38.83 

9%(la=3.5%) 

— 

— 


(1) Based upon analytic simulation and shake test data. 

(2) Computed based upon visual inspection of amplitude decay response 
resulting from removal of forced excitation. 

(3) Computed from shake test data. 

(4) Computed from Ct, “ 

K r 00^ 

(5) Forced excitation not used at 250 RPM because of large observed 
resonance at this condition. 

(6) Values shown are obtained at slightly different RPM. 
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Figure 1. 



Filter Desired 
Modal Response 


Identify Coefficients 
of Characteristic 
Equation 


Calculate Z-Transform 
Domain Roots 


Calculate Continuous 
Time Damping 
and Frequency 


0 )^, and 

Standard Deviations 


- Overall Procedure for Real Time Identification 
of Modal Damping and Frequency. 
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Figure 3- - Identified Parameter Convergence ( RJIS of Random Response = .25 deg, 
a^=.l deg, ^^“*2 deg/sec^, At=.02 sec, 0)^=23.15 rad/sec, ?=10%) . 


30 


Number of Samples 


Figure 4. 


Identified Parameter Convergence (RMS of Random Response = *05 
a =.l deg, a -.0^ deg/sec^, At=.02 sec, O) =23.15 rad/sec, ^-107o) 
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Figure 7. - Identified Parameter Convergence (RMS of Random Response 
a =.l deg, O =.2 deg/sec^, At=.02 sec, co =23.15 rad/sec, 
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Figure 11. - Identified Parameter Convergence Using Digital Filter 
Number 4 (RMS of Random Response = ,75 deg, ^^“-1 
O' =.2 deg/sec^, At=.02 sec, (o =23.15 rad/sec, ^=2%) * 
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Figure 12. - Effect of Filter Bandwidth on Identified Damping Ratio 
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Figure 15 


Hover Facility Test Data Used to Identify Damping and Frequency of the Inplane Mode at 
Q = 250 and = 618 RPM. (Only a Portion of the Data is shown, Amplitude Scale is Approximate) 
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Figure 17. - Identified Damping and Frequency Parameter Convergence From 
Hover Facility Test Data at fJ = 250 RPM (Parameter Accuracy 
is Best before Region of Bad Data is encountered at 1200 
samples, At=.02 sec) - 
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Figure 18. - Hover Facility Test Data at fi = 618 RPM Showii 
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Results are presented from simulation random response data v;hich quantify the 
identified parameter convergence behavior for various levels of random excitation. 

The data length required for acceptable parameter accuracy is shown to depend upon 
the amplitude of random response and the modal damping level* Random response 
amplitudes of 1.25 degrees to .05 degrees are investigated- The RML technique is 
applied to hingeless rotor test data from the NASA Langley Research Center Helicopter 
Hover Facility. The inplane lag regressing mode is identified at different rotor 
speeds. The identification from the test data is compared v/ith the simulation results 
and with other available estimates of frequency and damping. 
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